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Abstract 

The propagation of elastic waves in a fractured rock is investigated, both theo- 
retically and numerically. Outside the fractures, the propagation of compressional 
waves is described in the simple framework of one-dimensional linear elastodynam- 
ics. The focus here is on the interactions between the waves and fractures: for this 
purpose, the mechanical behavior of the fractures is modeled using nonlinear jump 
conditions deduced from the Bandis-Barton model classicaly used in geomechanics. 
Well-posedness of the initial-boundary value problem thus obtained is proved. Nu- 
merical modeling is performed by coupling a time-domain finite-difference scheme 
with an interface method accounting for the jump conditions. The numerical ex- 
periments show the effects of contact nonlinearities. The harmonics generated may 
provide a non-destructive means of evaluating the mechanical properties of frac- 
tures. 

Key words: elastic waves, contact nonlinearity, Bandis-Barton model, jump 
conditions, finite-difference schemes, interface method. 
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1 Introduction 



Fractures are the breaks in rocks caused by the huge stresses resulting from 
plate tectonics. It is of fundamental importance for geophysicists to be able to 
determine the position and the properties of fractures (such as their thickness) 
to be able to make predictions about the mechanical properties of a fractured 
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platform or the diffusion of a pollutant, for instance. Elastic waves are com- 
monly used for this purpose. When the wavelengths are much larger than the 
distance between fractures, the fractures are generally not studied individu- 
ally, and homogenization theories are applied. Otherwise, as in the case of 
the present study, it is possible to study single fractures. If, in addition, the 
wavelengths are much larger than the thickness of the fractures, the latter can 
be modeled in terms of interfaces with appropriate jump conditions. 



Many experimental, theoretical and numerical studies have dealt with wave 
propagation across thin and single fractures in terms of linear jump conditions 
[10,11,7]. The linear framework provides an appealing approach but it may not 
be very realistic, since non-physical penetration of both sides of the fractures 
may occur. Some authors have proposed more accurate fracture models, such 
as the Bandis-Barton model [2]. This model is commonly used in rock me- 
chanics and engineering to deal with quasistatic loading conditions, such as 
those occuring in flood barriers, for instance. However, very few studies have 
dealt so far with wave propagation in a model of this kind: we have only found 
one study using this appoach in [14]. 



The aim of the present paper is to further study the propagation of mechanical 
waves across a fracture described by the Bandis-Barton model, both theoret- 
ically and numerically, in a highly idealized configuration. To perform the 
numerical modeling, we combine an interface method with a classical numeri- 
cal scheme for wave propagation, namely the ADER scheme [13] . The interface 
method involves changing this scheme locally, using the jump conditions at 
the interface; moreover, it gives a subcell resolution when the interface does 
not coincide with the meshing. Many interface methods have been proposed 
since the 90's; see [6] for a review. Here, we adapt the explicit simplified in- 
terface method (ESIM) previously developed for dealing with linear contacts 
[7]; see [9] for an overview of the principles and advantages of this method. 
To our knowledge, this is the first time nonlinear jump conditions have been 
studied using an interface method 



The paper is organized as follows. In section 2, the problem under study is 
introduced, especially the hyperbolic jump conditions (4). An analysis of the 
solution is performed in section 3: conservation of energy, and the existence, 
uniqueness and regularity of the solution. Section 4 deals with the numerical 
methods. The numerical experiments performed in section 5 with realistic 
parameters show the influence of the incident wave amplitudes. 
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2 Problem statement 

2. 1 Configuration 

Consider a rock with a single plane fracture. Outside the fracture, the media 
involved f2j (i = 0, 1) are linearly elastic and isotropic; they are subject to a 
constant static stress —a (a > 0) running perpendicular to the fracture. At 
rest, the fractured zone is an interphase with thickness h > (figure 1, left). 
The nonlinear mechanical behavior of the interphase is investigated in section 




Fig. 1. Static (left) and dynamic (right) behavior of the fractured rock. I: incident 
wave; R: reflected wave; T: transmitted wave. 

Let us now consider a plane compressional wave propagating through Q nor- 
mal to the interphase; the interactions between this incident wave and the 
interphase give rise to reflected (in Q ) an d transmitted (in f^) plane com- 
pressional waves. These perturbations in fl and fli are described by the simple 
one-dimensional elastodynamic equations [1] 

dv da da 2 dv . . 

p ~dt = ~dx~ ) ~dt =pC 9^' [) 

where v = |f is the elastic velocity, u is the elastic displacement, and a is the 
elastic stress perturbation around —a. The physical parameters involved are 
the density p and the elastic speed of the compressional waves c; these piece- 
wise constant parameters may be discontinuous around the fracture: (p 0) 0)) if 
x G VIq, (pi, Ci) if x G Q\. The dynamic stresses induced by the elastic waves 
affect the thickness h(t) of the interphase (figure 1, right). Due to the finite 
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compressibility of the interphase, the constraint 

h = h + [u] > h - d > (2) 

must be satisfied, where [u] = u + — u~ is the difference between the elastic 
displacements on the two sides of the interphase, and d > is the maxi- 
mum allowable closure [2]. We also assume that the wavelengths of the elastic 
perturbations are much larger than h. One can therefore neglect the propaga- 
tion time through the interphase, and replace it by a zero-thickness interface 
at x — a, where a belongs to the interphase; therefore, [u] = [u(a, t)] = 
u(a + , t) — u(a~ , t). 

2.2 Bandis-Barton model 
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Fig. 2. Sketch of the stress-displacement relation deduced from the Bandis-Barton 
model. The bold straight segments denote the limit-case of unilateral contact. 

Single thin fractures have been classically modeled in terms of linear jump 
conditions [10]. Given a stiffness K > and neglecting the inertial effects, the 
most usual linear jump conditions are 

[a(a, t)} = 0, 
1 

[u(a, t)] = x a ( a ' 

The simple jump conditions (3) can be rigorously obtained by performing an 
asymptotic analysis of the wave propagation process within a plane interphase 
which is much thinner than the wavelength (/i < A); then K = pc 2 /h, where p 
and c are the physical parameters of the interphase. For K — > +oo, we obtain 
perfectly-bonded conditions; for K — > + , we obtain a(a ± , i) — > 0, and hence 
the two media VL Q and fli tend to be disconnected. The main drawback of the 
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conditions (3) is that they do not satisfy (2) under large compression loadings: 
a(a~, t) < —Kd =>- h < h — d, which contradicts (2). Hence, (3) is realistic 
only in the case of very small perturbations. When larger ones are involved, 
nonlinear jump conditions are required. 

To satisfy (2), we use the Bandis-Barton model [2]. This model is based on 
quasi-static compressional experiments showing that the closure of a fracture 
depends hyperbolically on the stress applied. In the case of dynamic problems 
[14], the hyperbolic jump conditions can be written 



with a(a~ , t) < K d (the other branch of the hyperbola is not realistic). In the 
Bandis-Barton model [2], the parameters K and d are linked to a, satisfying 
a < K d. The second relation in (4) is sketched in figure 2. Under compression 
loadings: a(a ± , t) < 0, the second equation of (4) implies: [u(a, t)] > —d, 
hence (2) is satisfied. Under traction loadings: <r(a ± , t) > 0, (2) is trivially 
satisfied. In the latter case, nothing prevents a(a~, t) > a, leading to dis- 
connection or adhesion processes; more realistic models that account for such 
processes are not investigated here. The straight line with a slope K tangen- 
tial to the hyperbola at the origin describes the linear jump conditions (3); 
as deduced from (4), the linear conditions are valid only if |it(q; ± , t)\ <C K d. 
In the limit-case d — > + with K bounded, the hyperbola tends towards the 
nondifferentiable graph of the unilateral contact, denoted by bold straight 
segments in figure 2. The latter graph corresponds to the Signorini's condi- 
tions [12]: <r(a + , t) = a(a~, t) < 0, [u(a, t)] > 0, and cr(a ± , t) [u(a, t)] = 0. 
This limit-case, which leads to a difficult mathematical analysis and requires 
suitable numerical tools, is not investigated here. 



2.3 Initial-boundary value problem 

From now, we follow a velocity-stress formulation of elastodynamics. It re- 
mains therefore to know the jump condition satisfied by v at a. For that 
purpose, we differentiate the second equation of (4) with respect to t 



[a{a, t)} 







[u(a, t)] 



1 a(a , t) 
K a(a~, t) 
Kd 



(4) 
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and we replace the time derivative in (5) by a spatial derivative via (1). To- 
gether with the first equation of (4), it gives 



v (a + , t) = v(a , t) + 



cr(a + , t) = a (a , t). 



K 



dv 



1 - 



a(a~, t)\ d 



x 



a , t), 



Kd 



(6) 



To sum up the conservation laws (1) together with the jump conditions (6), 
we set up 

/ 

A(x) 



U(x, t) 



W 







p 



\-pc 



2 



(7) 



and then we state the following initial-boundary value problem 

x 7^ a, t >t , 



( d d 
— U + A — U = forxe 
at ox 



U(a + , t) = D (u(a-, t), ^U(a~, t)^ , 
U(x, t ) = Uq(x) for x eR, 



(8) 



where Dn : M 4 



is a nonlinear application deduced from (6) and (7). We 



assume that the initial data Uq(x) : 
support included in fl , with p > 1. 



>2 ; 



is a function with a compact 



For use in section 4.2, we need also the jump conditions satisfied by some 
of the spatial derivatives of U (under suitable assumptions of regularity as 
defined in theorem 1). Differentiating (6) m — 1 times with respect to time, 
and then replacing the time derivatives by spatial derivatives via (1), yield 
nonlinear m-th order jump conditions that can be written 

U(a + , t)=D m [ U{oT, t), w-^U(a-, t), ——U{oT, t) , (9) 



d x 



dx r ' 



d x 



where D m : R 2(m + 2 ) — > R 2 is a nonlinear application. The computation of 
D rn is tedious task that can be automated using computer algebra tools. The 
simulations shown in section 5 were obtained in this way. 



3 Analysis of the solution 



The aim of this subsection is to prove that the initial-boundary value problem 
(8) is a well-posed problem. Our proof is based on rather elementary tools, 
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and yields an analytical solution. The numerical evaluation of the latter with 
a fourth-order Runge-Kutta scheme constitutes the semi-analytical solution 
used in section 5. Before giving the theorem, we need some intermediate re- 
sults. 



3.1 Conservation of energy 

In the first lemma, we define an energy, and we show that it is conserved. 
Lemma 1 Let U(x, t) be a solution of (8). Then, 



PROOF. We multiply the first equation of (1) by v, and we integrate it by 
parts on fl ; then the second equation of (1) gives 




(10) 



satisfies 




0, E(U, t) > 0, E(U, t) = U(x, t) = 0. 



dt 





d 




In the same way, we obtain 
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Adding the two previous equations gives: S\ + e 2 = 0, where 



£1 



r+00 Q v r+oc I Q a 

/ p v — — ax + — 2 a dx 
J — 00 '<) I .1 



pc 2 dt 



d 1 f +oc ( 2 1 2 ^ 
= — - / pv H ^ a ax, 



Ex 



and, using (5), 



£2 



a(a , £) [f (a, £)] , 

1 <7(aT, t) 9(7 
^7 t7(a-, t)\ 2 9l 



(a , t), 



dt 



( 



In 1 - 



cr(a! , t) 
Kd 



a (a , £) 
Kd 



- 1 



£9 



The quantity E — Ei + E 2 therefore satisfies = 0; -Ei is obviously a 
positive definite quadratic form. All we have to do now is to study the sign of 
E 2 . Setting 9 = l-a(a~, t)/(Kd), a study of g(9) = hi 9 + 1/9 - 1 for 9 > 
shows that g{9) > 0, and that g{9) = = 1, i.e. a{pT , t) = 0. From (4) 
and (12), one sees that <r(ct~ , t) = o-(q; + , t) and u (a 1 * 1 , t) = 0. □ 

Lemma 1 means that (10) is an energy which is split into two terms. The first 
term, Ei, is the classical mechanical energy associated with the propagation 
of elastic waves outside the fracture. The second term, E 2 , is the mechanical 
energy associated with the nonlinear deformation of the fracture. Since the 
inertial effects are neglected in (4), E 2 amounts to a potential energy. Note 
that in the limit case \<r(a~, t)\ <C K d, one gets E 2 — ■> 2 -^<t 2 (q;~, t), which 
corresponds to the well-known potential energy of a linear spring. 



3. 2 Method of characteristics 



To express U(x, t) in terms of limit- values of the fields at a, we now use the 
Riemann invariants J R,L that are constant along the characteristics 7^ [5]. 
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The invariants for linear PDE's with constant coefficients are very simple 



1R ■ 


dx 

Tt =+C = 


dJ R 
> — — 
dt 


= 0, 


with J R (x, t) 








1L ■ 


dx 

— = -c = 
dt 


d.J L 
dt 


= o, 


with J L (x, t) 



1 v -yA^^ 



1L 



v H a (x, t). 



(11) 

After some calculations, the first equation of (4) and (11) along with the initial 
data condition of compact support in Q give for t > t 



a(a ± , t) = —p\ c\v(a + , t), 



v(a , t) 



Pi Cl 



:i2) 



v(a + , t) + 2 J* (a - c (t - t ), t ) , 



Poc 

T R,L 



where the subscript i on J i ' refers to The following lemma expresses 
U(x, t) in terms of v(a + , s), with to < s < t, and of the initial values of the 
Riemann invariants, which are linked to the initial data U (x) via (7), (8) and 
(11). 

Lemma 2 Setting 

1 / 1 / 
tA — t [a — x), ts —t (x — a), 

C Ci 



the solution U(x, t) of (8) is given by 
( , , \ / 



x < a : U(x, t) 



1 1 

y -Po c Po c J 



Jo(x-c (t-to), t ) 
A A (x, t) 



with Aa(x, t) 



' 11 v(a + , t A ) + J R (a - c (t A - t ), t ) if t A > t , 



Po c 



x > a : U(x, t) = 



Jq(x + co(t — t ), t ) otherwise, 

( i \ 



-Pi Ci 



A B (x, t), 



with Ab(x, t) = < 



v(a + , t B ) if t B > t , 
otherwise. 



(13) 
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Proof. The properties of Riemann invariants (11) in Qq lead to 



Jq{x, t) — Jq 



f(x - co(t - t ), h), 
Jo{a~, t A ) if t A > t , 
Jq(x + c (t — t ), t ) otherwise. 



(x, t) = < 



For (x, t) such that t A > t , (11) and (12) imply that 

Jo(a-, t A ) = -P±^ v ( a + ,t A ) + J^a - c (t A - t ), t ). 
Po c 

The value of and the previous equations allow to conclude for x < a. Since 
the support of U (x) is included in f2 , (H) in ^i lead to 



otherwise, 
Jf(x, t) = JL(x + a(t- t ), t ) = 0. 

For (x, t) such that t B > t , (11) and (12) imply that 

(a + , t B ) = v(a + , t B ). 

The value of A B and the previous equations allow to conclude for x > a. □ 

Lemma 3 The limit value y = v(a + , t) satisfies the nonlinear ODE 



J*(x, t) = 



Jf (a+ t B ) if t B > to, 



-jj = f(y, t), y(to) = o, with 



^'j-^o^'H^-^sS')- (14) 



g(t) = 2J R (a-co(t-t ), t ) . 



PROOF. From (6) and the first equation of (12), we deduce 




1 



1 



da 



Pi Ci 
K 



1 



dv 



(a + , t). 




v(ct , t) is then eliminated via the second equation of (12), giving (14). 



□ 
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3.3 Well-posedness of the initial-boundary value problem 



As shown by lemmas 2 and 3, the solution of (8) is expressed in a unique 
manner in terms of v(a + , t) solution of the ODE (14). Proving the existence 
and uniqueness of the solution to (8) therefore amounts to showing that the 
solution to (14) exists and is unique, as followed. 

Theorem 1 Let Uq{x) G C p (R), p > 1. There exists a unique global solution 
U(x, t) G C p (R \ a, [to, +oo[) to the initial-boundary value problem (8). 

PROOF. For t > t ,y — > f(y, t) in (14) is C°°, it is therefore a locally Lipschitz 
function. Moreover, Jq(x, t ) is a C v c function (p > 1) deduced from U (x), 
therefore g(t) in (14) is C v c : hence t — > f(y, t) is continuous. The Cauchy- 
Lipschitz theorem ensures that the solution y(t) is unique, if it exists. Lastly, 
/ G C p implies that y G C p+1 [4]: the lemma 2 ensures the C p regularity of U. 

Suppose that y(t) is not bounded as t — > t*; the first equation of (12) im- 
plies that a(a~, t) — > ±oo. Since a(a~, t) < K d (see section 2.2), only the 
case o~(a~, t) — > — oo needs to be addressed. In this case, (10) implies that 
E{U , t) — > +oo, which is impossible: lemma 1 and U G C p mean that 
E(U, t) = E(U , t ) < +oo. Hence y(t) is always bounded, and the local 
existence due to the Cauchy-Peano theorem is also global [3]. □ 



4 Numerical treatment 



4-1 Numerical scheme 





a 


a 


+ 




X 


J 


-i j 


J 


+1 J 


— 5^ 

+2 



Fig. 3. ID rock fractured at x = a; spatial mesh. 

Given (xj, t n ) = (i Ax, to + n At), where Arc is the mesh size and At is 
the time step, we seek an approximation U™ of U(xi, t n ). We use two-step, 
explicit, and (2 s + l)-point spatially-centered finite-difference schemes to in- 
tegrate (8). Time-stepping can then be written symbolically 



Ur 1 = U? + H q (UI S ,...,UQ, 



(15) 



with q = if Xi G tt , q = 1 otherwise, and with H q : R 2x ( 2s+1 ) -> R 2 [5]. 
We define J so that xj < a < xj + i (figure 3). A grid point is regular if all 
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the grid points Xi- S , Xi +S used in (15) belong to the same medium as x«; 
in this case, (15) is applied classically. Otherwise, a grid point is irregular, 
and its time-stepping is described in the next section. The irregular points are 

x .J-s+li ■■■ 1 X.J+S- 

For the numerical experiments in section 5, we choose the ADER r schemes 
[13], where r denotes the order of accuracy. The stability limit is CFL = 
max(c) At / Ai = 1. For odd values of r, as chosen here, s = r / 2. Other 
schemes can be used (numerical experiments have been performed successfully 
with a flux-limiter scheme and a fifth-order WENO scheme): a priori, readers 
can adapt their favorite solver to the forthcoming discussion. 



4-2 Interface method 



One applies the explicit simplified interface method (ESIM) [7,9]: at irregular 
points, some of the numerical values used for the time-stepping procedure (15) 
are modified. At time t n and at irregular points Xi, these modified values U* 
are numerical estimates of smooth extension U*(x, t n ) of U(x, t n ) across a 



(16) 



2k ~ 1 (x - a) m d m 
x>a, U*(x,t n )=Y,— -j ^ U ( a 

m=0 ■ UX 

2k-l I _ \m am 

x < a, E/*(x, t n ) = £ 1 " J ' , m U(a + , t n ). 

^ ml ox m 

The next paragraph details how to estimate J^c U(a ± , t n ) in (16). 

Estimation of U(a~, t n ). We focus here on the computation of the mod- 
ified values at irregular points Xj on Q\ (j = J+l, J+s); the case of irregular 
points on f2 (j — J—s + 1, J) is similar. To estimate U(a~, t n ) in (16), 
we write Taylor expansions of U (xj, t n ) on the left of a (i = J — k + 1, J) 

2fc-l / _ \m am 

U(x it t n )= £ ^U(a-,t n ) + 0(Ax 2k ), (17) 

m =o • ° x 

and on the right of a (i — J + 1, J + with the jump conditions (9) 

2fc-l C r . _ 3m 

C(*«*») = E ^^J^(«Vn) + 0(Ax 2fc ) 



m=0 m ! <9x r 

2fe— 1 / .Am / am 



1 x,- — a 



Qra Qm+1 



m=0 



+0(Ax 2fc ). 

(18) 
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The exact values U(xi, t n ) are replaced by the known in (17) and (18); 
■^p-k U(a~,t n ) and the Taylor remainders are also removed, leading to 

Flu',-, U-) = (lT}_ k+1 , U% k ) T , (19) 



where -Jj^ U are estimates of -J^ U(a , t n ), and where i* 1 : M 4fc — >■ M 4fc is a 
nonlinear application. The system (19) is solved using Newton's method. 

Computation and use of modified values. Once the limit values -Jj^ U~ 
are known, the modified values U* are deduced from U*(x i: t n ) (16): 



2fc-l 

E 

m=0 



i=j + u...,j + ., ^;=E^r. (20) 



When all the modified values have been computed, all that remains to be done 
is to perform time-stepping at irregular points. To do so, the time-stepping 
(15) at an irregular point X{ requires the use of the numerical values at the 
grid points in the same medium as Xi, and the modified values otherwise 

i = J - s + 1, .., J, = U? + H (Ut s , U n j, U* J+1 , U* +s ) , 

i = J + l,...,J + s, C/™ +1 = 17? + H 1 (U*_ s , U*j, U n J+1 , U^ +s ) . 

(21) 

Comments. We do not propose a numerical analysis of the interface method 
in the nonlinear context. We just recall some properties that are true in the 
limit case of negligible nonlinearities; readers are referred to [7] for proofs. For 
\<r(a~, t)\/(K d) = 0, the following properies are satisfied: 

(i) the system (19) has always a unique solution; 

(ii) in the limit case of a homogeneous medium without any fracture (i.e., 
Po Pi, c o —> c i) K +oo) and if k > s, then U* — > J7": the interface 
method (21) amounts to the classical time-marching (15); 

(in) for a r-th order scheme, the local truncation error of (21) is still r-th 
order accurate if 2 k — 1 > r. 

The convergence measurements performed in section 5 indicate that the prop- 
erty (in) is still satisfied in the nonlinear context, but a rigorous proof remains 
to be obtained. Up to now, we have not chosen k: for that purpose, we fol- 
low two criteria. First, the spatial derivatives in (17) and (18) need to be 
well-defined: theorem 1 implies 2 k — 1 < p. Second, we want the properties 
(i,ii,iii) to be satisfied in the linear limit. It leads to the following inequalities 



13 



2 k - 1 < p, 
< k>s, (22) 
2k-l>r. 

In practice, we use the minimum value of k that satisfies (22). 

We have no theoretical results about the stability of the interface method, 
even in the linear context. We have studied many geometrical configurations, 
values of the physical parameters, and values of K, d. With a wide range of 
parameters, no instabilities are usually observed up to the CFL limit, even 
after very long integration times. However, instabilities are observed when the 
physical parameters differ considerably between the two sides of an interface. 
In practice, this is not penalizing when dealing with realistic media. 

The nonlinear system (19) may have more than one solution, and we can not 
be sure that Newton's algorithm selects the right one. To check this point, 
we compare the numerical value v~ obtained by solving (19) with the semi- 
analytical value deduced from (14) and (12). For weak to moderate nonlinear 
effects, the numerical value and the semi-analytical value are the same (up 
to the accuracy of the integrations). But when the nonlinear effects are large, 
convergence towards a wrong solution can occur if the mesh used is too coarse. 
This is not very surprising: since wave profiles tend to be stiffened (see section 
5), the Taylor expansions in (17) and (18) give rather poor estimates. The 
strategy followed in section 5 consists in using a finer mesh. A better strategy 
might consist in solving (19) with the constraint a~ < K d, as required by (4). 



5 Numerical experiments 

We study a 400-m domain fractured at a — 200.67 m, with parameters [14] 
Po = Pl = 1200 kg/m 3 , K = 1.3. 10 9 kg/s 2 , 

< 

c = Cl = 2800 m/s, d = 6.1 10~ 4 m. 

Since p and c are the same on both sides of a, the reflected wave and the 
distorsion of the transmitted wave are induced only by the mechanical behavior 
of the fracture. Many numerical experiments have been tested successfully with 
different physical parameters on both sides of a. 
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5.1 Test 1: initial-boundary value problem 

In the first experiment, the initial data is U (x) = (— l/c , po) T h(t — x / c ), 
where h is a spatially-bounded C® combination of sinusoids 

UE 8 -" sin(/3 m u; c £) if < f < — , 
I otherwise, 

with (3 m = 2 m ~ 1 , oo c = 2ir f c ; the coefficients a m are: ai = 1, a 2 = —21/32, a 3 = 
63/768, a 4 = —1/512. The central frequency is f c = uj c /(2n) = 50 Hz, and 
to = 52 ms. This choice ensures that the incident wave is a purely rightward- 
travelling wave, originally located in Qq. Three values of e are considered, 
leading to three amplitudes t>o of v(x, to): 0.01 m/s (a), 0.1 m/s (b), and 0.4 
m/s (c). ADER 4 is coupled with the ESIM with k = 3, as defined by (22). 
The computations are performed with CFL = 0.9 on N = 400 grid points ((a) 
and (b)) or N = 1200 grid points (c): here the number of grid points is larger, 
otherwise (19) would not have converged to the right solution. 

In the left column of Figure 4, one shows the numerical and analytical values 
of a at t — 116.29 ms. At this instant, the incident wave has crossed the 
fracture and one sees the reflected and transmitted waves. In the right column 
of Figure 4, one shows the time history of the numerical values of [u(a, t)] 
deduced from (4) and (19). Here the horizontal dotted line represents —d, and 
the vertical scales are the same for (a), (b) and (c). 

In case (a), Vq is too small to mobilize the nonlinearity of the fracture: almost 
no differences could be detected with simulations performed with (3). Case 
(b) corresponds to realistic seismic waves recorded during on-site investiga- 
tions: moderate nonlinear effects are present. Case (c) corresponds to incident 
blasting waves: large nonlinear effects are present, stiffening the fronts. In the 
right column, one clearly sees how min([w]) approaches — d when the incident 
wave amplitudes increases, without reaching this value as required by (2). 

5.2 Test 2: convergence measurements 

To check the accuracy of the the interface method, we compare the analytical 
and numerical values of a on successive refined meshes. The parameters are 
the same than in subsection 5.1. The errors are measured in norm L 1 at t — 
116.29 ms. All the computing is carried out in double precision. 
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Fig. 4. Test 1: v = 0.01 m/s (a), v = 0.1 m/s (b), and v = 0.4 m/s (c). Left 
column: numerical (...) and exact (-) values of a; right column: numerical values of 
[u] (the horizontal dotted line denotes —d). 
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Fig. 5. Test 2: convergence measurements for ADER 2 (a) and ADER 4 (b). 

Figure 5 shows convergence measurements with ADER 2 (a) and ADER 4 (b) 
coupled with the ESIM ((22) implies k = 2 for ADER 2). The values of v 
are the same as in test 1: 0.01 m/s (o), 0.1 m/s (*), and 0.4 m/s (•). The 
vertical scale of (b) is twice as large than the vertical scale of (a). In all cases, 
the expected orders of accuracy are exactly reached: order 2 in the case of 
ADER 2, order 4 in that of ADER 4. These results indicate that coupling a 
r-th order accurate scheme with the interface method is still r-th order if the 
appropriate value of k (22) is used. 



5.3 Test 3: sinusoidal source term 



As a last test, we simulate a time-harmonic experiment. To do so, we slightly 
change the system under study: instead of (8), we solve 

f d d 

— U + A — U = 5(x-x s ) S(t) for x ^ a, t>0, 

' U(a\ t) = D (u(a~, t), ^U(a~, f)) , ( 24 ) 

U(x, 0) = for x G K, 

where S(t) = (0, e sm{oj c t)) T is a source of mass that acts at x s = 40 m (in 
Q ). Except for the source, the parameters are the same than in subsection 5.1. 
Three values of e are considered, yielding the three amplitudes v : 0.01 m/s 
(a), 0.1 m/s (b), and 0.4 m/s (c). The periodic values of the transmitted wave 
are recorded at x — 220 m; a decomposition into Fourier series is then applied 
to these values. Figure 6 shows the numerical values of a on one period (left 
column) and the normalized coefficients of the Fourier series (right column). 
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Fig. 6. Test 3: vo = 0.01 m/s (a), vq = 0.1 m/s (b), and vq = 0.4 m/s (c). Left column: 
numerical values of a during a period; right column: normalized coefficients of the 
Fourier decomposition. 
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The harmonics generated when vq increases are clearly seen in this figure. 

There are two reasons for displaying these harmonics. First, they help to decide 
whether it is worthwhile taking the nonlinear effects into account: solving the 
nonlinear system (19) at each time step is more costly than solving a linear 
system during a pre-processing step, as done with the linear jump conditions 
(3) in [7]. In case (a), the answer is negative; in case (b), it depends on the 
accuracy required; in case (c), the answer is positive. Secondly, the decrease 
in the harmonics is linked to the stiffness K and the maximum allowable 
closure d. It may therefore be possible to infer these values by inspecting the 
harmonics, which is a non-destructive means of evaluating the fracture. 



6 Conclusion 

Here we have studied the propagation of 1-D elastic compressional waves across 
a contact nonlinearity. The latter feature models a fracture in rocks, but it 
can also be a useful means of describing other physical situations, such as 
those encountered in the nondestructive evaluation of material for instance. 
The present study involves a physical description of the model, a theoretical 
analysis of its solution, and a numerical time-domain modeling. 

The numerical analysis of the interface method needs to be studied further to 
determine its accuracy and stability (the latter point is still an open question 
in the linear context). The resolution of the nonlinear system (19) needs also 
to be improved. Lastly, it would be interesting to determine the harmonics 
generated across the fracture analytically, using the harmonic balance method. 

Many difficulties need to be overcome before tackling realistic 2-D configura- 
tions. First, it is required to model the mechanical shear behavior of realistic 
fractures, along with the compressional behavior [2]. These effects can be effi- 
ciently modeled in 1-D configurations. Secondly, the 2-D treatment proposed 
in [8] to deal with linear contacts is greatly intricated: the linear underdeter- 
mined systems of jump conditions are now nonlinear. Thirdly, the computa- 
tional extra-cost induced by the interface method is substantially higher: the 
extrapolation matrices used in the interface method need to be computed at 
each time step and at many grid points along the fracture. Fourthly, if the 
nonlinear effects are large, a fine grid will be required to converge towards the 
right solution of the nonlinear systems. Using a local mesh refinement around 
the fracture might therefore be a useful strategy here. 
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